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ABSTRACT 

For  the  band-limited  inverse  problem,  the  inversion  formulas  developed  by 
Bleistein  et  ah  generate  a  reflector  map  as  well  as  an  estimate  of  reflection 
coefficients  from  the  discontinuities  of  the  medium  parameters.  Here,  we  de¬ 
scribe  the  use  of  a  multiscale  smoothing  operator  to  suppress  the  noise  while 
preserving  the  features  of  the  original  processing  formula.  This  is  accomplished 
by  applying  a  simple  multiplier  in  wavenumber  domain  to  the  output  of  the 
inversion  formalism.  A  comparison  of  this  approach  to  wavelet  based  edge  de¬ 
tection  processing  shows  that  this  continuous  multiscale  operator  allows  more 
flexibility  for  changing  the  length  of  the  smoothing  operator. 

Key  words:  multiscale  smoothing  operator,  noise  suppression,  Fourier  inver¬ 
sion,  imaging,  wavelet  transform 


Introduction 

Rapid  changes  in  the  medium  parameters  characterize 
the  structure  of  the  subsurface;  they  are  the  reflectors  in 
the  earth.  In  situations  where  the  medium  parameter  is 
related  to  a  velocity  field,  or  the  perturbation  in  the  ve¬ 
locity  field,  it  is  typical  to  think  in  terms  of  piecewise 
smooth  functions  whose  discontinuity  surfaces  represent 
reflectors.  The  object  velocity  fields  discussed  in  this  re¬ 
port  all  belong  to  this  class  of  step-like  piecewise  smooth 
functions. 

The  inversion  formulas  developed  by  Bleistein  et 
0^.(1996)  provide  a  tool  for  obtaining  correct  locations 
of  interfaces  cis  well  as  model-consistent  specular  reflec¬ 
tion  coefficients.  However,  because  of  the  band-limited 
nature  of  seismic  data,  the  step-like  wavespeed  perturb¬ 
ation  functions  are  not  well  reconstructed  using  this  for¬ 
mulation.  By  modifying  the  inversion  operator  (Bleistein 
et  ah,  1996),  the  output  can  be  transformed  into  the  sin¬ 
gular  functions  of  the  discontinuity  surfaces  of  the  ori¬ 
ginal  function,  scaled  by  the  jump  in  the  original  function 
at  each  point  of  the  discontinuity  surface.  This  modifica¬ 
tion  is  a  normal  derivative  operator  in  the  spatial  domain, 
achieved  through  an  appropriate  multiplier  in  the  Fourier 
domain  derived  from  asymptotic  analysis. 

The  derivative  operation  is  sensitive  to  the  noise 
present  in  the  wavefield.  In  this  report,  a  multiscale 


smoothing  procedure  is  introduced  to  stabilize  the  inver¬ 
sion  algorithm  and  improve  the  images.  The  multiscale 
operator  is  based  on  the  Fourier  transform.  The  al¬ 
gorithm  extracts  the  rapid  changes  in  the  velocity  field 
as  a  function  of  location  and  scale,  while  noise  is  largely 
suppressed. 

More  specifically,  we  introduce  a  scale,  s,  in  the  con¬ 
volution  operator.  It  is  known  (Mallat  &:  Hwang,  1992), 
that  the  effect  of  the  scaling  operator  on  noise,  Ar(x),  is 
to  produce  an  output  Ns  (x)  having  the  property  that 

|iV3(x)|  =  ©(s"),  u<0.  (1) 

That  is,  Ns{x.)  decays  with  increasing  s.  On  the  other 
hand,  for  a  piecewise  smooth  signal,  /(x),  we  know  that 
the  scaled  outputs  satisfy 

|J,(x)|  =  0{s°).  (2) 

Thus,  its  peak  amplitude  does  not  decay  with  increas¬ 
ing  s.  However,  when  /(x)  is  a  delta  function,  /s(x)  is 
a  bandlimited  delta  function  whose  resolution  decreases 
with  increasing  s.  Thus,  it  is  a  “race”  between  noise  sup¬ 
pression  and  resolution.  Our  tests  with  parameters  typ¬ 
ical  of  seismic  data  suggest  that  we  win  the  race  in  this 
application. 

The  Bleistein- Cohen  (BC)  inversion  formulas  have 
the  structure  of  Fourier  transform-like  integrals.  One  fun¬ 
damental  concept  of  the  BC  inversion  approach  is  ap- 
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plying  the  underlying  relationship  between  the  observa¬ 
tions  and  the  Fourier  transform  of  the  function  (s)  that 
characterize (s)  the  unknown  medium.  Thus,  in  the  fol¬ 
lowing  sections,  attention  will  be  primarily  focused  on 
a  multiscale  smoothing  operator  applied  to  Fourier  ima¬ 
ging  as  a  test  of  its  ability  to  suppress  noise.  Thereafter, 
we  address  the  applications  of  this  method  to  the  BC 
inversion. 

This  discussion  adds  the  structure  of  scaling  analysis 
to  the  fairly  standard  noise  reduction  processes  currently 
in  use. 


Fourier  Imaging 

In  seismic  applications  such  as  filtering,  imaging,  and  mi¬ 
gration/inversion,  one  basic  concept  is  the  reconstruction 
of  the  data  by  performing  forward  and  inverse  Fourier 
transforms  with  proper  “focusing”  filters  that  depend  on 
both  domains.  We  axe  concerned  with  aperture-limited 
Fourier-like  transforms  of  data  /(x').  Such  an  output  can 
be  written  as  a  cascade  of  forward  and  inverse  Fourier 
transforms  acting  on  /(x'), 

J(x)  =  ^  frk  f  dV  (3) 

V  )  J Dj^  J D^t 

where  n  is  the  dimension  of  interest;  x  and  x'  are  n  com¬ 
ponent  vectors;  k  is  the  wave  vector.  D^f  and  Dk  are  the 
supports  of  /(x')  in  space  and  wavenumber  domain  re¬ 
spectively;  that  is,  they  define  the  apertures.  In  seismic 
data  processing,  Dk  is  constrained  to  a  “large  wavenum¬ 
ber”  domain;  the  nature  of  the  Dk  is  essential  to  the  re¬ 
construction  of  /(x').  In  the  applications  to  BC  inversion, 
/  is  related  to  the  wavespeed  perturbation  of  the  sub¬ 
surface.  It  represents  a  piecewise  smooth  function,  from 
which  it  follows  that  the  domain  of  the  integration  can  be 
decomposed  into  separate  domains  whose  boundaries  in¬ 
clude  all  of  the  discontinuities  of  /.  The  amplitude  of  the 
integrand  a(x',k)  is  the  “focusing”  operator  required  to 
reconstruct  the  image  from  the  data;  it  is  allowed  to  de¬ 
pend  on  both  k  and  x'.  The  appropriate  filter  arises  from 
properties  of  aperture-limited  large  wavenumber  Fourier 
transforms.  To  simplify  the  problem  and  focus  our  atten¬ 
tion  on  some  insights  of  Fourier  inversion,  we  start  our 
discussion  with  a  “  1  in  (3). 

Bleistein  et  a/.  (1996),  have  shown  that  the  high 
wavenumber  aperture-limited  Fourier  inversion  of  the 
piecewise  smooth  function  /(x^)  is  dominated  by  the  val¬ 
ues  on  the  discontinuity  surfaces.  The  inversion  is  ap¬ 
proximately  a  band-limited  step  function  of  normal  dis¬ 
tance  from  each  point  on  the  discontinuity  surface.  The 
amplitude  of  the  step  function  is  proportional  to  the  jump 
in  the  function  across  the  surface  at  the  point  in  question. 


For  better  imaging  of  the  discontinuity  surfaces, 
the  following  filter  in  wavenumber  domain  is  applied 
(Bleistein  et  al,  1996): 

zA;sgn(u  •  k).  (4) 

That  is,  the  filter  a(x',  k)  is  replaced  by  the  above  multi¬ 
plier.  Here,  k  =  |k|  is  the  magnitude  of  vector  k,  k  is  the 
unit  vector  and  u  is  a  constant  vector.  In  seismic  applic¬ 
ations,  there  is  usually  little  or  no  information  about  the 
plane,  kz  =  0,  so  u  can  be  chosen  to  be  (0,  ±1)  in  2D  or 
(0, 0,  ±1)  in  3D,  respectively,  so  that  sgn(u  k)  =  dhsgnfcz. 
The  multiplication  by  iksgn{u-k)  replaces  the  Fourier  in¬ 
version  of  /(x')  by  the  singular  function  (s)  of  its  bound¬ 
ary  surface (s),  with  the  amplitude  proportional  to  the 
jump  in  /(x')  at  each  point  on  the  surface.  The  method 
breaks  down  at  points  where  u  •  n  is  nearly  zero,  with  n 
being  the  unit  normal  to  the  reflector.  In  the  application 
to  inversion,  iksgn{u  ■  k)  is  replaced  by  iuj]  v  ^|)  with 
(j)  being  traveltime.  In  this  case,  no  direction  u  is  dis¬ 
tinguished;  all  directions  are  treated  the  same.  Bleistein 
et  a/.  (1996)  have  shown  that  this  filter  behaves  asymp¬ 
totically  like  a  normal  derivative  operator  in  the  spatial 
domain,  even  though  the  normal  direction  is  not  known 
a  priori. 

Multiscale  Operator  for  Noise  Suppression 

This  directional  derivative  operator  is  sensitive  to  noise. 
It  enhances  noise  in  the  data  at  leirge  wavenumbers  more 
than  at  small  wavenumbers.  In  order  to  extract  the  im¬ 
portant  object  boundaries  from  the  data,  we  need  to 
find  the  local  maxima  representing  the  peaks  of  the 
singular  functions  while  minimizing  the  effect  of  noise. 
A  multiscale  smoothing  operation  can  do  the  trick.  It 
extracts  the  most  rapid  changes  in  the  neighborhood 
defined  by  the  scale  parameter  while  largely  suppressing 
the  noise. 

The  smoothing  process  can  be  viewed  as  a  convolu¬ 
tion  operation, 

e(x')  +  /(x').  (5) 

Define  0s  (x')  as  the  multiscale  transform,  given  by, 

©s(x')  =  0(x7s)-  (6) 

Here,  s  is  the  scale  parameter.  It  is  used  to  control  the 
size  of  the  neighborhood  where  the  local  maxima  of  deriv¬ 
atives  are  computed.  At  coarse  scales  -  large  s  -  a  large 
neighborhood  is  encompassed,  wherecis  at  fine  scales  - 
small  s  -  small  neighborhoods  contribute  to  the  convolu¬ 
tion. 

In  this  report,  we  choose  the  smoothing  operator 
0(x')  to  be  a  Gaussian,  i.e.. 
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Figure  1.  (a)  Velocity  model,  (b)  velocity  singularities,  (c)  velocity  singularities  of  noisy  data  with  signal  to  noise  ratio  of  4/3, 
and  (d)  image  of  (c)  by  applying  the  smoothing  operator. 


©(x')  =  exp  (— x'  ■  x  /2).  (7) 

The  Gaussian  filtering  kernel  function  is  infinitely  differ¬ 
entiable,  and  has  the  same  form  in  wavenumber  domain. 
Thus,  the  operation  is  just  the  following  multiplication, 

0s(k)  =  V^sexp  (— s^k  •  k/2),  (8) 

in  the  wavenumber  domain  before  inverting  the  trans¬ 
form. 

By  combining  the  multiscale  smoothing  operator 
with  the  directional  derivative  operator  in  the  previous 
section,  we  obtain  the  multiscale  operator, 

^^(k)  =  zA;sgn(u  •  k)05(k).  (9) 


Replacing  a(x',  k)  in  (3)  with  (9),  we  propose  the  follow¬ 
ing  processing  formula  for  detection  of  the  discontinuities 
of  /(x'),  while  simultaneously  suppressing  the  noise: 

=  /crfc4-.(k)  f  dV 

V  JDk 

Figure  1  is  an  example  of  applying  this  multiscale 
singularity  detection  procedure  in  two  dimensions.  Fig¬ 
ure  1-a  shows  a  2D  earth  model  with  constant  velocity 
in  layers  and  strong  discontinuities.  Rapid  changes  in  the 
medium  parameters  are  clearly  indicated  in  Figure  1-b. 
Here,  we  have  applied  the  normal  derivative  operator  (4) 
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Figure  2.  (a)  Single  trace  of  full-band  data,  (b)  band-limited 
data,  (c)  singularities  of  full-band  data  (a),  (d)  singularities 
of  band-limited  dada  (b),  (e)  band-limited  noisy  data  and  (f) 
image  of  (e). 


to  the  noise-free  data.  When  noise  is  added  to  the  original 
model  at  a  signal  to  noise  ratio  of  4/3,  the  directional  de¬ 
rivative  operator  (4)  enhances  the  noise  as  is  shown  in 
Figure  1-c.  The  structural  information  is  seemingly  lost, 
however,  the  rapid  changes  can  still  be  recovered  at  a 


certain  scale  s,  using  (8-10).  This  is  shown  in  Figure  1-d 
where  we  have  computed  (9)  with  s  =  32.  The  scaling 
behavior  of  the  noise  is  such  that  it  is  largely  suppressed 
when  compared  with  Figure  1-c.  This  result  compares 
favorably  with  Figure  1-b. 

There  are  various  algorithms  based  on  wavelet  trans¬ 
forms  that  address  this  goal  of  noise  suppression  (Song 
&  He,  1996),  (Dessing  et  aL,  1996).  Basically,  all  the 
algorithms  consist  of  an  operator  having  a  behavior 
roughly  comparable  to  a  (band-limited)  differentiation 
and  smoothing  operation  followed  by  extraction  of  the 
local  maxima  of  the  derivative  of  the  data.  The  local  max¬ 
ima  occurs  at  the  singular  surfaces  and  in  the  direction 
perpendicular  to  them.  The  direction  of  the  preference  is 
not  known  a  priori.  Thus,  the  wavelet  operation  is  often 
carried  out  by  computation  of  the  modulus  function  from 
mutually  perpendicular  wavelet  transforms.  The  modu¬ 
lus  function  is  essentially  the  magnitude  of  the  gradi¬ 
ent  operator.  Thus,  it  can  extract  the  local  maximum 
changes  within  the  region  defined  by  the  scale  parameter; 
however,  it  fails  to  indicate  the  sign  of  the  data  changes. 
On  the  other  hand,  the  multiplier  (4)  will  produce  a  nor¬ 
mal  derivative  except  when  that  normal  is  itself  perpen¬ 
dicular  to  u*.  Because  seismic  reflectors  are  rarely  ab¬ 
solutely  vertical,  we  can  choose  u  as  a  unit  vector  in 
the  vertical  direction  and  image  the  normal  derivative  in 
a  practical  range  of  normal  directions.  An  advantage  of 
this  operator  is  that  it  reveals  the  magnitude  of  the  step 
jumps  by  its  spike  amplitude.  The  white  sinusoid  in  Fig¬ 
ure  1(d)  shows  that  it  has  a  negative  amplitude,  which 
is  consistent  with  the  decrease  of  the  wavespeed  in  the 
original  model. 

For  further  illustration,  we  process  one  trace  of  the 
output  in  Figure  1.  Figure  2-a  shows  the  trace.  The 
Nyquist  wavenumber  is  50/A;m.  Figure  2-b  shows  the  de¬ 
gradation  of  the  bandlimited  2  —  16/ km  data  when  com¬ 
pared  to  the  full-band  data  (Figure  2-a).  The  effect  of  los¬ 
ing  low  wavenumber  components  is  the  flattening  of  the 
data  away  from  the  discontinuities;  if  the  data  at  zero  fre¬ 
quency  is  zero,  then  the  average  of  the  inverse  wavenum¬ 
bers  must  be  zero.  Further  low  wavenumber  filtering 
transforms  the  steps  into  doublets.  The  high  wavenum¬ 
ber  filtering  only  affects  the  width  of  the  doublets;  16 /km 
causes  insignificant  loss  of  resolution.  The  steps  are  un¬ 
determinable  although  the  locations  of  the  discontinuities 
of  the  data  are  recognizable.  Figure  2-c  shows  the  singu¬ 
larities  of  the  full-band  noise-free  data  (Figure  2-a);  it 
is  the  output  of  applying  the  same  procedure  as  shown 


*  However,  by  the  observation  in  the  first  section,  that,  in  ap¬ 
plication  in  our  inversion  technique,  there  is  no  distinguished 
direction,  u. 
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in  Figure  1-b.  Figure  2-d  is  the  corresponding  result  for 
band-limited  data  2-b.  A  comparison  of  Figure  2-c  and 
2-d  implies  that  the  bandlimiting  effect  on  the  normal  de¬ 
rivative  operator  (4)  is  only  the  insignificant  differences 
on  the  widths  of  resolutions  determined  by  the  band- 
widths.  The  degradation  of  bandlimiting  on  the  original 
data  has  been  overcomed  by  employing  the  normal  deriv¬ 
ative  multiplier  (4).  Figure  2-e  is  the  band-limited  noisy 
data  by  adding  noise  to  the  band-limited  data  2-b  at  a  sig¬ 
nal  to  noise  ratio  of  4/3.  Figure  2-f  shows  the  result  of  ap¬ 
plying  the  noise  filtering  derivative  operator  proposed  in 
the  paper  to  these  data.  By  comparing  the  peak  locations 
here  with  those  in  Figure  2-b,  we  see  that  the  discontinu¬ 
ities  are  still  well  recovered  from  the  band-limited  data, 
even  though  the  low  wavenumber  components  that  con¬ 
tain  information  about  the  steps  of  the  original  function 
are  missing.  The  peaks  of  the  Gaussians  clearly  reveal  the 
locations  of  the  discontinuities,  while  their  heights  indic¬ 
ate  the  magnitude  of  the  step  scaled  by  the  area  under 
the  filter  applied  to  the  data.  The  widths  of  the  Gaussi¬ 
ans  are  a  manifestation  of  the  high  pass  filtering  effect 
of  setting  s  =  32;  this  choice  of  scale  was  determined  as 
a  compromise  between  higher  values  through  which  fur¬ 
ther  resolution  was  lost  and  lower  values  for  which  the 
noise  level  was  not  acceptable. 


Singularity  Reconstruction  and  Amplitude 
Preservation 

In  this  section,  we  will  analyze  some  features  of  the 
multiscale  Gaussian  operator  (9)  as  applied  to  singularity 
detection.  We  will  show  that  this  operator  reconstructs 
the  discontinuities  of  the  function  by  Gaussians.  The  loc¬ 
ation  of  each  discontinuity  is  indicated  by  the  peak  of  the 
Gaussian,  while  the  step  height  proportional  to  the  peak 
amplitude.  The  amplitude  remains  unchanged  across  dif¬ 
ferent  values  of  the  scale  s. 

We  consider  one  specific  piecewise  smooth  function 
of  /(x')  in  (10):  /(x)  =  AH{v  •  (x  —  xo)),  where  R(x) 
is  a  step  function  having  discontinuities  at  the  set  of  xo; 
V  denotes  the  normal  at  each  point  on  the  discontinuous 
surface;  and  A  is  a  scalar.  When  in  (10)  includes  all 
the  discontinuities,  and  Dk  is  of  infinite  extent,  the  calcu¬ 
lation  can  be  carried  out  in  a  similar  way  as  in  the  corres¬ 
ponding  aperture-limited  Fourier  inversion  examples  in 
(Bleistein  et  aL,  1996).  Here  we  simply  state  the  solution, 

/,(x)  =  (11) 

This  implies  that  the  full-band  multiscale  Fourier  inver¬ 
sion  characterizes  the  singularities  by  the  Gaussian.  The 
peak  of  Gaussian,  when  x  =  xo,  indicates  the  location  of 


the  singularity,  and  its  peak  amplitude  is  just  the  mag¬ 
nitude  of  the  step  jump.  It  will  not  change  with  varying 


Let  us  now  consider  the  effect  of  bandlimiting  on 
these  results.  We  assume  that  Dk  is  a  region  symmetric 
about  k  =  0,  so  that  we  can  then  consider  symmetric 
intervals  in  each  kj^j  —  l,2,...7i.  The  solution  can  be 
rewritten  as. 


Is(x) 


exp{- 


(X  -  Xq)]^  ' 
2s2  ^ 


+T]{S,X-  Xo,Kt,_,Ki;+  )]  ,  (12) 

where  Kv  is  the  wavenumber  in  the  v  direction;  Kv_  and 
are  the  two  end  points  of  7y(s,  x  —  xq,  ,  «v+) 
is  a  lower  order  term  and  equals  0  at  x  =  xo  where  the 
maximum  of  |/s(x)|  occurs. 

Equation  (12)  indicates  the  effect  of  band-limited  in¬ 
version  compared  to  the  full-band  Fourier  image.  The 
difference  between  (11)  and  (12)  at  their  peak  values  is  a 
difference  of  error  functions  governed  by  the  limits  of  in¬ 
tegration  in  Kir.  For  appropriate  choices  of  kv±  and  s,  the 
bandlimited  output  in  (12)  is  indistinguishable  from  the 
full  bandwidth  version  in  (11).  Again,  the  above  equation 
implies  that  the  peak  amplitude  is  the  step  jump  scaled 
by  a  factor  that  is  asymptotically  constant  in  s,  namely 
1. 


The  Scale  Parameter  and  Characterization  of 
the  Singular  Behavior  of  Functions 

The  behavior  of  the  multiscale  operator  (9)  is  such  that 
it  characterizes  the  velocity  singularities  within  the  re¬ 
gion  defined  by  the  scale  parameter  while  smoothing  the 
region.  In  mathematics,  singularities  are  generally  char¬ 
acterized  by  their  Lipschitz  exponents.  A  function  f{x) 
has  Lipschitz  exponent  of  i/  at  xq  if  and  only  if  there  ex¬ 
ists  a  constant  B  such  that  for  all  x  in  a  neighborhood 
of  Xo,  we  have 

l/(®)  -  /(s^o)l  <  B\x  -  xol",  v>v.  (13) 

That  is,  is  a  lower  bound  on  the  powers  of  differences 
in  X  that  bound  differences  in  the  function  values.  The 
function,  /(x),  is  uniformly  Lipschitz  v  over  an  open  in¬ 
terval  (a,  6)  if  and  only  if  the  constant  B  holds  for  any 
X,  Xo  €  (a,  h).  For  example,  if  /(x)  is  a  step  function  dis¬ 
continuous  at  Xo,  then  its  Lipschitz  exponent  is  z/  =  0. 
The  isolated  singularity  is  the  worst  singularity  inside 
a  region  that  contains  the  isolated  point.  The  uniform 
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Lipschitz  regularity  of  a  function  over  a  certain  region 
is  then  equal  to  the  pointwise  Lipschitz  regularity  at  the 
isolated  point.  (Song  Sz  He,  1996). 

When  noise  is  present,  the  observed  data  becomes 
/(x)  H-n(x),  and  the  image  ^(x)  in  (10)  becomes  /^(x); 
that  is, 

I,(x)  =  /3(x)  +  iV,(x).  (14) 

If  n(x)  is  white  noise,  the  Lipschitz  exponent  i/  of  Ns  (x) 
is  less  than  zero  (Song  &  He,  1996),  (Mallat  &  Hwang, 
1992).  That  is,  there  exists  a  constant  such  that 

\Ns{^)\<Bs^,  u<0,  (15) 

This  inequality  implies  that  Ns  (x)  should  decrease  when 
the  scale  s  increases.  On  the  other  hand,  the  peak  val¬ 
ues  of  Is  remains  unchanged  with  s  increasing  according 
to  equations  (11)  and  (12).  Hence,  depending  upon  the 
differences  between  the  signal  singularities  and  the  noise 
singularities,  we  can  take  the  advantage  of  the  multiscale 
operator  (9)  to  suppress  the  noise  in  the  data  at  a  faster 
rate  than  the  loss  of  resolution  of  the  derivative  at  dis¬ 
continuities.  We  also  remark  that  the  scale  s  should  stay 
bounded  above  for  seismic  data,  because  of  the  spatial 
resolution  being  lost  for  increasing  scale. 

The  scale  parameter  determines  the  size  of  the  win¬ 
dow  for  optimal  peak  detection  combined  with  noise  sup¬ 
pression.  In  different  problems,  we  can  consider  changing 
the  window  by  simply  varying  the  scale  parameter  within 
a  reasonable  range.  In  the  pure  discrete  wavelet  approach, 
scaling  is  usually  discretized  to  powers  of  2  to  facilitate 
synthesis  of  the  processed  signal  from  its  wavelet  trans¬ 
form.  In  contrast  the  multiscale  operator  has  the  advant¬ 
age  that  we  search  s  values  over  a  continuous  range  to 
optimize  our  choice  for  simultaneous  singal  detection  and 
noise  suppression.  In  this  application  of  scaling  analysis, 
there  is  not  need  for  synthesis  ~  we  will  extract  all  in¬ 
formation  from  the  transform  of  the  original  data.  Hence, 
we  need  not  restrict  our  scales  to  powers  of  2. 

Applications  to  BC  Inversion 

The  multiscale  smoothing  operator  is  essentially  a 
smooth  low-pass  filter  in  wavenumber  domain.  It  facilit¬ 
ates  the  numerical  identification  of  the  discontinuity  sur- 
face(s)  and  the  amplitude  of  the  jump  at  each  point  on  the 
surface(s),  while  suppressing  the  noise  in  a  large  scale. 
Thus,  we  can  apply  this  procedure  to  Bleistein’s  inver¬ 
sion  formula, 

q(x)  =  [ 

•  J  dcjexp  {— iu;(^(x,^)}us(x5,  Xsjo;).  (16) 


Here,  Ck(x)  is  the  perturbation  correction  to  the  given 
reference  velocity  field;  the  spatial  weighting  B(x,f)  con¬ 
tains  a  determinant  that  characterizes  the  viability  of  in¬ 
verting  source-receiver  configurations;  Us  is  assumed  to 
be  the  observed  data  from  a  medium  characterized  by 
piecewise-smooth  parameter  functions  with  discontinu¬ 
ity  surfaces  behaving  as  reflectors.  With  the  band-limited 
high-frequency  seismic  input  data,  the  output  of  equation 
(16)  is  a  band-limited  step-like  function,  with  its  discon¬ 
tinuity  and  amplitude  not  very  well  recovered.  The  noise 
in  the  observed  data  is  propagated  into  a(x)  through 
the  procedure.  The  multiscale  operator  can  then  be  ap¬ 
plied  as  a  posterior  process.  The  output  is  a  series  of 
band-limited  delta  functions.  The  locations  and  the  re¬ 
flection  coefficients  can  be  extracted,  while  the  noise  is 
suppressed  because  of  the  properties  of  this  multiscale 
operator,  as  shown  above. 

Conclusions 

We  have  described  here  a  multiscale  smoothing  operator 
that  simultaneously  suppresses  noise  and  accurately  de¬ 
tects  the  location  of  discontinuities  of  a  piecewise  smooth 
function  of  the  type  that  represents  velocity  models  in 
seismic  inversion.  Analytical  and  numerical  investiga¬ 
tions  confirm  these  results  and  predict  that  the  size  of  the 
discontinuity  can  be  estimated  from  the  output.  These 
properties  are  exactly  what  is  needed  in  the  implement¬ 
ation  of  this  operator  in  the  BC  inversion  formalism,  but 
now,  with  the  added  feature  of  noise  suppression. 
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